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PREFACE 


The  research  discussed  in  this  report  had  severe!  objectives:  1) 

an  examination  of  the  effect  of  the  free  surface  on  the  near-field  seismic 
radiation  for  seismic  sources  imbedded  in  an  elastic  half-space;  2)  an 
examination  of  the  interrelation  of  various  parameters  such  as  source  dimen- 
sion, rise  time  and  frequency  and  3)  an  examination  of  the  effects  of  fault 
length  and  rupture  velocity  in  the  near-field  together  with  a study  of 
the  near-field  coda. 
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FIGURE  10. 
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Assumed  fault  model  and  location  of  observing  stations. 
Transverse  component  of  displacement  U^  at  Stations  A and  D. 
Amplitude  is  given  in  terms  of  the  assumed  U^  in  cm  divided 
by  1000.  T is  the  rise  time  in  seconds  and  v is  the  rupture 
velocity. 

Transverse  component  of  displacement  U^  at  stations  A1  and 
Dl.  A1  and  D1  are  located  at  distances  less  than  the  critical 
distance  of  the  SP  head  wave. 

Three  components  of  displacement  at  Station  B located  in  the 
SW  quadrant  in  the  direction  of  fault  rupture. 

Three  components  of  displacement  at  station  C located  in  the 
NW  quadrant  behind  the  direction  of  fault  rupture. 

Displacement  seismograms  at  stations  B and  C which  illustrate 
the  effect  of  varying  the  rupture  velocity  for  fixed  source 
rise  times  of  0.5  and  1.0  seconds. 

Displacement  seismograms  at  stations  B1  and  Cl  for  various 
rupture  parameters. 

Vertical  component  of  displacement  at  stations  B and  C for 
various  rupture  parameters. 

Vertical  component  of  particle  velocity  at  station  B for 
various  rupture  parameters.  Note  the  increase  in  amplitude 
as  the  rupture  velocity  approaches  the  shear  velocity  (3.2 
km/sec)  of  the  medium. 

Vertical  component  of  particle  acceleration  at  station  B for 
a source  rise  time  of  0.5  seconds.  Note  the  appearance  of 
the  P and  S wave  stopping  phases. 
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FIGURE  11. 


Vertical  component  of  particle  acceleration  at  station  B for 


FIGURE  12. 


FIGURE  13. 


FIGURE  14. 


FIGURE  15. 


FIGURE  16. 


FIGURE  17. 


a source  rise  time  of  1.0  seconds.  Wave  character  is  simpler 
than  in  Figure  10  and  individual  stopping  phases  are  not 
recognizable  since  the  source  time  is  comparable  to  the 
fault  length  divided  by  the  rupture  velocity. 

Vertical  component  of  particle  acceleration  at  station  C.  Note 

» 

the  distinct  stopping  phases. 

Vertical  component  of  displacement  for  a buried  strike-slip  source 
at  two  epicentral  distances. 

Radial  component  of  displacement  for  a buried  strike-slip  source 
at  increasing  epicentral  distances.  Note  appearance  of  Rayleigh 
pulse . 

Radial  component  of  displacement  for  a buried  dip-slip  source 
at  increasing  epicentral  distances.  The  seismogram  for  r = 100  km 
begins  at  16.5  seconds  and  has  been  offset  to  avoid  overlapping 
seismogram  at  r = 50  km.  Amplitude  scales  are  identical. 

Radial  component  of  displacement  for  a dip-slip  source  at 
increasing  epicentral  distances.  Note  enhancement  of  Rayleigh 
pulse  due  to  shallower  depth  of  source  as  compared  to  Figure  15. 
Seismogram  at  r = 5 km  utilizes  displacement  scale  ranging  from 
0 to  -20  microns.  Seismograms  at  10,  20,  50  and  100  km  are 
plotted  at  same  amplitude  scale  utilizing  scale  next  to  r = 10  km. 
Seismogram  at  r = 100  km  begins  at  16.5  seconds  and  slightly 
overlaps  end  of  r = 50  km  seismogram. 

Observed  and  computed  near-field  seismograms  for  a central 
California  earthquake.  A is  the  fault  area  and  u is  the  average 
dislocation.  Different  values  of  Au  are  estimated  on  opposite 
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sides  of  the  fault  presumably  due  to  the  effect  of  geology 


at  the  recording  sites.  SAGO-East  is  situated  on  Pliocene 
sediments  whereas  SAGO-Central  is  located  on  granite.  Larger 
amplitudes  are  commonly  observed  on  sediments  relative  to 
bedrock. 

FIGURE  18.  Theoretical  point  source  seismograms. 

a)  Homogeneous  half-space. 

b)  Layered  half-space.  Source  depth  12.5  km,  epicentral 
distance  5 km. 
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INTRODUCTION 
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Observations  have  recently  been  made  with  broadband  displacement 
transducers  recording  within  a few  source  depths,  source  dimensions  or 
wavelengths  of  central  California  earthquakes  (Johnson  and  McEvilly,  1974). 

At  these  short  distances,  where  hopefully  the  effects  of  scattering  and 
attenuation  are  minimized,  one  hopes  to  gain  insight  into  details  of  the 
faulting  mechanism,  such  as  rupture  velocity  and  the  time  function  of 
displacement . 

Our  motivation  for  computing  theoretical  seismograms  is  to  gain 
some  insight  into  the  effects  of  the  free  surface,  and  rupture  parameters 
for  a propagating  strike-slip  fault  in  the  near-field.  Specification  of  a 
working  model  and  its  resultant  theoretical  seismogram,  which  can  be  com- 
pared to  observed  data,  is  an  important  facet  in  our  ability  to  predict 
strong  ground  motion  in  the  near  field  of  potentially  damaging  earthquakes. 

Near-field  displacements  have  been  examined  using  dislocation  modeling 
in  an  infinite  space  (Haskell,  1969;  Kanamori,  1972;  Trifunac,  1974;  Trifunac 
and  Udwadia,  1974;  Anderson,  1974;  Anderson  and  Richards,  1975).  The  effect 
of  a free  surface  on  near-field  displacements  has  been  examined  by  Kawaski 
et  al.  (1973)  and  Kawasaki  and  Suzuki  (1974).  Anderson  (1976)  has  examined 
the  near-field  motions  due  to  a shallow  rupturing  fault  using  a Green's 
function  which  is  a solution  to  Lamb's  problem.  Anderson's  paper  criti- 
cally discusses  the  conditions  under  which  the  free  surface  may  be  reason- 
ably accounted  for  by  doubling  the  amplitude  of  motions  in  an  infinite  space. 
In  particular,  the  SH  motions  from  a strike-slip  source  and  the  P-SV  motions 
from  a dip-slip  source  were  examined.  Archuleta  and  Frazier  (1976)  have 
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utilized  a three-dimensional  finite  element  approach  to  examine  particle 
displacements  and  velocities  in  the  near  field  for  propagating  shear  frac- 
tures. 

In  this  research,  we  have  made  use  of  the  generalized  multipolar  ray 
tiieory  (GMRT)  advocated  by  Ben-Menahem  and  Vered  (1973)  to  examine  the 
near-field  motions  of  a propagating  strike-slip  fault.  GMRT  is  a comprehen- 
sive theory  adaptable  to  the  computation  of  both  near-field  and  far-field 
theoretical  seismograms  utilizing  surface  or  buried  dislocation  sources. 

The  theory  also  permits,  under  certain  conditions,  the  computation  of  ground 
velocity  and  acceleration. 


THEORY 


We  adopt  a cartesian  coordinate  system  in  which  the  free  surface  of  a 


homogeneous  half  space  is  represented  by  the  xy  plane.  In  the  xz 


plane  we  define  a rectangular  fault  surface 


0 < X < L 


Hf  < z < 


Displacement  occurs  such  that  all  points  of  the  rectangle  having  the 


same  coordinate  x move  simultaneously,  i.e.  a line  source  with  a width 


h = H.,  - H moves  along  the  xz  plane. 


It  has  been  shown  by  Ben-Menahem  and  Vered  (1973)  that  the  Laplace 


transformed  surface  displacements  of  each  comnonent  due  to  a dislocation 


point  source  located  at  (0,0, H)  are  given  by  expressions  of  the  form 


— 1 ;Cos  .M(fi  n ’ , , 2£+nr(-l^,  2.  , / 2 1 ^1/2,.,,  ... 

U ^ (sur)u  f(u  )exp[-s(u  + ) H]du  (1) 


N sin  M(f) 


£,  n,  m,  N and  M are  source  dependent  integer  constants.  r is  the 


epicentral  distance,  ({)  is  the  azimuthal  angle,  c is  the  elastic  wave 


propagation  velocity,  s is  the  Laplace  parameter,  u is  an  inte- 


2 2 

gration  variable  and  f(u“')  is  a function  of  u . A step-like  source 


time  function  is  initially  assumed  at  the  source,  i.e..  '-'q(s)  = U^^s 


The  depth  of  the  source,  H,  appears  only  in  the  exponent  in  (P. 


Thus,  using  the  principle  of  superposition,  the  Laplace-transformed  dis- 


placement generated  by  a line  source  which  is  compi'sed  o\  point  sources 


continuously  distributed  along  an  interval  h = H , - H , is  given  b'  tin. 


express  ton 


1 n-1  fCOS  M(^ , , , ^ 2£-Hn+l  f(u  ) 

UdH  = -^  s { . i J (sur)u  

, N sin  M4>  ■ m'  " 1/2 

1 o'  . ^ . -L 


(u  + -2) 


•{exp[-s(u^  + - exp[-s(u^  + )^^^H2])du 

c c 


(2) 


We  denote  by  U^,  and  the  vertical,  radial  and  azimuthal 


components  of  displacement.  From  Table  2 in  Ben-Menahem  and  Vered 
(1973)  and  (2)  we  obtain 


- ^ dx  sin  24i  ^ 

z 4tt 


r 


( — E - aE,  )du 

" " D 


0 f A ^ a a 

j 


(3) 


- ^ dx  cos^_  y 


27T 


r 

' < I + 72  "b>  - ^ ■*" 

b b 


(4) 


2F 


^r  = zd^LJinJi  ^ ( T E,  + -|  E J + (b ' E,  - ^ E J } du  (5) 


4TT 


srA  ' a a ' 2 b"  A 

b 


a a 


where 


-saHj^  -saH2 
[e  - e 


-sbH  -sbH„ 

f i.  ^ ' 

[e  - e : 


m+1 

u J (sur) 
m 


P = fb’  - 2ab] 


A = b'^  - u^ab 
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o - / 2 ^ , 2, 1/2 

a - (u  + 1/a  ) 


b = (u^  + 1/3^)^^^ 


b'  = (u^  + 1/23^) 


The  longitudinal  and  shear  wave  velocities  are  denoted  by  a and  6 
reflectively. 

The  time  domain  solution  is  obtained  by  the  method  of  Cagniard- 
Pekeris  for  dislocation  sources.  As  an  example  we  will  treat  the  radial 
component  of  motion.  We  first  write  (5)  in  the  form 


U =U  +-U  +U  +-U 
r r.pj^  s r.p^  r.s^^  s r,S2 


(6) 


where 


dx  sin  2({) 


r u\ 

— E du 
aA  a 


(7) 


r,p. 


-dx  sin  2( 
2irr 


, r b 

-4-  E du 
aA  a 


(8) 


-dx  sin  2({) 
Ait 


Fb’ 


(9) 


r ,s. 


^ r p 

-dx  sin  2 <t>  ■ 2 

“TT— * "o  -J-,  h 

. b A 


(10) 
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where 


IL 

2 ^ 2 
R C 


2-1/2  R 

T ) T < — 

' ’ c 


I£+  iH  (,2  _ ^ 1/2  ^ R 

2 2^  2 ' ’ r 

R R C 


2 2 2 
R = r + H 


C = a or 


. r R r ^ , 1 1 ,1/2- 

min[  e , a + H ( -^  - — ) ] 

6 a 


T = R/B 
s 


= f(H^)  - f/H^) 


From  (6)  we  thus  obtain 


U^(t)  = U (t)  + U (x)dT  + U (t)  + ' U 

^ r,p^  r,s^  j r.J 


(t  )dT 


Equation  (15)  gives  the  ground  displacement  due  to  a step  function  U 


For  a continuous  source  function  U^gft)  we  denote  the  ground  displacement 
by  G^(t).  G^(t)  is  given  by 
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t 


G^(t)  = U^(t-T)g'  ( r)dr 


G^(t)  = U|^(T)g(t-T)dx 


From  (15)-(17)  we  obtain 


G (t)  = {[U  ^ (T)  + U ^ (x)]g’(t-T)  + 


[U  (x)  + U (x) ]g(t-x) } dx 


By  assuming  that  g(t)  has  first  and  second  continuous  derivatives  we 
can  obtain  ground  velocity  and  acceleration  by  differentiating  (18): 


c;(t) 


{ [U  (x)  + U (x)  ]g"  (t-x)  + 


[U  (t)  + U (x)  ]g' (t-x)  Idx 


G';(t)  = 


U (x)  + U (x)g’  " (t-x) 


[U  (x)  + U (x)]g"(t-x)l  dx 


Table  1 gives  the  necessary  change  of  variables  to  evaluate  the 
integrals  in  (11)-(14),  Calculated  seismograms  for  the  infinite  media  can 
be  accomplished  by  a slight  change  of  the  integrands  in  (11)-(14). 
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Table  1 


Integral  Transformations 


q 


sin 


r 

c 


^62-’ 


q 


. -1 
Sin 


r 

c 


q = sin 


1/2 


c 

r 

s 


T 

S 


write 


T 


S 


q 


. -1,1  a/2 

sin  ( — ) 


To  obtain  seismograms  for  a finite  source,  we  next  divide  the 


r 


I 


length  of  the  fault  L into  equal  subintervals  6L^  (i  = 1,  r) . 

For  every  subinterval,  the  corresponding  time  function  g(t)  derived 
from  expression  (1)  or  (2)  is  computed  for  the  selected  azimuthal  angle 
t and  epicentral  distance  r.  The  final  signal  is  calculated  by 
adding  the  signals  corresponding  to  all  the  subintervals  taking 

into  account  the  appropriate  time  delay.  A source  time  function 


g(t) 


r 0 


t < 0 , 

0 ^ t T,  (21) 

t > T , 


where  t is  the  adopted  rise  time  is  assumed.  The  total  time  of 
rupture  is  equal  to  (L/v)  + t where  L is  the  fault  length  and  v 
is  the  rupture  velocity.  The  assumed  source  time  function  is  basically 
a step  function  with  rounded  shoulders  and  it  will  be  useful  for  purposes 
of  discussion  to  define  an  effective  rise  time  equal  to  0.6  t . 

The  effective  rise  time  is  simply  the  rise  time  of  a ramp  step 

function  positioned  along  the  time  axis  to  approximate  the  response  of 
the  rounded  step  function. 

The  numerical  results  were  checked  by  independently  computing 

the  static  displacements  for  a buried  strike-slip  point  source.  In  an 

elastic  half-space  the  static  displacements  are  given  by  the  expressions 

(Ben-Menahem  and  Singh,  1968): 

___  nisin2((i 

U = A tan^  — 


^ f!isin2ij>  tp 

U = [AW  - B]  cosecijj  tan^  — 

r2  2 


(22) 


14 


cos2(j)  4i 

[4a  - 2]  coseci);  tan^  — 

r2  2 


where  fii  = Uq  dS/Air, 

Uq  = amount  of  dislocation, 

dS  = fault  area, 

ij/  = angle  between  vertical  axis  through  source  and  radius  vector 
to  point  of  observation, 

W = cosij; 

a = Poisson's  ratio, 

A = [3(1  + W)2w  - (1  - 2a)(2  + W) ] , 

B = [3(1  + W)2  - |1  - 2a) ] . 

The  static  displacements  in  an  infinite  space  can  be  derived  by  making 
use  of  the  Laplace  transform  relation 

lim  U(t)  = lim  s U(s) 

t cx)  s -*•  0 

3S2i  sin2<j)  ^ 

U = cosij/  tan^(il;/2)  (1  + cosij;)^  , 

^ 4(l-a)R2 

__  fiisin2(j) 

U = (3  cos^ij/  + 4a  - 5)  cosecii'  tan^(4)/2)(l  + cosi|i)^  (23) 

^ 4(l-a)R2 

cos2([) 

= (2a  - 1)  cosecij;  tan^  (iii/2)  (1  + cosy)^. 

2(1  - a)R2 

All  displacement  components  for  all  sources  tend  to  a limit,  the  so-called 
"residual  deformation."  In  the  near-field  this  limit  is  reached  very 
soon  after  the  arrival  of  the  S-wave.  Evaluation  of  the  ratio  of  the 
half-space  static  displacements  to  those  of  the  infinite  space  for 
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various  angles  points  out  that  the  static  displacements  for  an  elastic 
half-space  cannot  in  all  cases  be  allowed  for  by  simply  doubling  the 
Infinite  space  displacements.  This  conclusion  is  in  agreement  with 
the  results  of  Anderson  (1976) . A further  check  on  the  infinite  space 
calculations  was  made  using  results  obtained  using  the  computer  program 
developed  by  Boatright  and  Boore  (1975). 


NUMERICAL  RESULTS 


The  coordinate  system  for  the  model  assumed  in  this  paper  is  shown  in 

Figure  1.  A P-wave  velocity  a = 5.5  km/sec  and  an  S-wave  velocity  6 = 3.2km/sec 

2 

has  been  used  throughout.  A rectangular  fault  surface  having  an  area  of  3 km 
is  imbedded  in  an  elastic  half-space  at  a depth  of  5 km.  5 km  can  be  taken 
as  a representative  depth  for  California  earthquakes.  The  fault  plane  is 
arbitrarily  taken  to  lie  in  a N-S  direction  with  rupture  initiating  at  the 
northern  end  and  propagating  in  a southerly  direction  axis) . We  shall 

be  examining  the  U^,  U2  and  U^  components  of  the  derived  displacements, 
velocities  and  accelerations  at  various  station  locations.  Stations  A1  a^  A 
lie  along  the  strike  of  the  fault  in  the  direction  of  rupture  and  stations 
B1  and  B at  an  angle  of  45°  to  the  strike  (azimuth  235°).  Stations  Cl  and  C 
are  located  at  an  azimuth  of  325°  and  stations  D1  and  D are  along  the  strike 
but  behind Jthe  direction  of  rupture.  In  the  results  which  follow  seismograms 
are  computed  for  rupture  velocities  of  2,  2.5  and  3 km/sec  and  source  rise 
times  T of  0.5  and  1 seconds.  The  computed  seismograms  are  given  in  units 
of  Uq/1000  cm  where  Uq  is  the  assumed  displacement  on  the  fault. 

We  first  proceed  to  examine  the  effect  of  the  inclusion  of  a free  surface. 

In  Figure  2 we  show  theoretical  seismograms  for  the  transverse  component  of 
motion  (U^)  at  station  A and  D with  and  without  the  presence  of  the  free  sur- 
face. As  can  be  seen  the  pulse  shapes  are  similar  and  the  peak  to  peak 
amplitudes  are  doubled  when  the  free  surface  is  included.  This  agrees  well 
with  theory  for  an  Incident  SH  component  of  motion.  Figure  3 compares  the 
transverse  component  of  displacement  at  stations  A1  and  D1  located  closer  to 
the  ends  of  the  fault.  At  these  locations  the  peak  to  peak  amplitudes  are 
only  approximately  doubled,  ranging  from  a factor  of  1.7  to  2,  when  the  free 
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surface  is  included.  Stations  A1  and  D1  are  located  at  distances  less  than 
the  critical  distance  for  the  existence  of  the  SP  head  wave. 

In  Figures  4 and  5 we  examine  the  effect  of  the  presence  of  a free 
surface  on  the  three  components  of  displacement  at  stations  B and  C.  As 
before  the  overall  pulse  shapes  are  quite  similar  and  the  main  effect  of 
tlie  presence  of  the  free  surface  is  to  increase  the  amplitude  of  the  wave 
pulse.  As  before  a factor  of  2 increase  in  amplitude  is  a good  approximation 
although  in  some  instances,  particularly  for  the  vertical  component  of  dis- 
placement , the  difference  in  amplitudes  is  closer  to  a factor  of  2.5.  It 
thus  appears,  that  for  the  strike-slip  example  presented  here,  allowance  for 
a free  surface  by  doubling  the  amplitude  of  the  wave  pulse  computed  for  an 
infinite  space  is  valid  to  within  ±25%  or  so.  Anderson  (1976)  has  pointed 
out  that  the  SV  contribution  to  the  displacement  is  strongly  influenced 
by  the  presence  of  the  free  surface,  depending  on  the  angle  of  incidence,  and 
is  particularly  emphasized  with  a dip-slip  source. 

As  can  be  seen  in  Figures  2-5  the  static  displacements  computed  with 
the  presence  of  a free  surface  differ  by  a factor  of  1 to  approximately  2 times 
the  infinite  space  static  displacements  depending  on  the  component  of 
motion  and  the  station  location.  It  can  be  aemonstrated  using  the 
expressions  given  by  (22-23)  that  this  is  indeed  the  expected  result,  i.e., 
the  static  displacements  cannot  in  all  cases  be  allowed  for  by  simply 
doubling  the  infinite  space  displacements. 

Figures  6,  7 and  8 show  displacement  seismograms  computed  at  stations 
B,  C,  B1  and  Cl  which  compare  the  effects  of  varying  the  rupture  velocity  for 
fixed  source  rise  times  of  0.5  and  1.0  seconds.  The  effect  of  a free  surface 
has  been  included.  For  a fixed  source  rise  time  increasing  the  rupture 
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velocity  from  2 to  3 km/sec  produces  a small  increase  in  amplitude  and  a 
slight  decrease  in  the  width  of  the  wave  pulse.  This  is  to  be  expected  inas- 
much as  the  width  of  the  pulse  should  be  a function  of  the  total  time  interval 
in  which  the  fault  is  rupturing.  If  we  define  as  the  distance  from  the 

edge  in  which  rupture  begins  to  the  observing  point  the  P-wave  initial  motion 

begins  at  R /a.  The  radiation  from  the  source  is  completed  at  a time  — + t ; 

a V e 

S-waves  initiating  at  this  moment  will  arrive  at  the  point  of  observation  at 
a time  where  R^^  is  the  direct  distance  from  the  edge  in  which 

rupture  is  completed.  The  time 


T 


(24) 


seems  to  control  the  pulse  width  of  the  near-field  seismograms  computed  here. 

ySf 

The  theoretical  seismograms  presented  here  all  possess  pulse  widths  equal  to  T. 
Equation  24  describes  the  dependence  of  rise  time,  rupture  velocity  and  the 
direction  of  rupture.  The  effect  of  the  direction  of  rupture  may  be  demon- 
strated by  comparing  seismograms  at  stations  A and  D and  A1  and  Dl.  For  ex- 
ample, the  difference  in  the  value  of  T for  stations  D and  A is  about  1 second 
and  for  stations  Dl  and  A1  is  about  0.6  second.  This  difference  can  be  seen 
on  the  computed  seismograms.  The  reason  that  the  pulse  width  is  simply  T 
appears  to  result  from  the  fact  that  the  static  displacement  approximately 
achieves  its  final  value  soon  after  the  arrival  of  the  S-wave  contributed 
from  the  segment  in  which  the  fracture  is  terminated.  This  will  in  general, 
however  not  always  be  true;  for  example,  at  larger  distances  where  the  emerging 
Rayleigh  wave  will  contribute  to  a larger  pulse  width  than  that  given  by  (24)  . 

For  a point  source  the  pulse  width  is  equal  to  ^ ^ ~ a ^ ' 
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It  may  In?  noted  that  the  terms  L/v  and  in  (24)  are  independent  of 

the  azimuth  to  the  station.  Thus,  one  may  expect  that  different  kinematical 

fault  models  possessing  the  same  L/v  + will  result  in  the  same  pulse 

width.  This  implies  that  rupture  velocity  is  not  separable  by  examining  the 

pulse  width  at  stations  of  varying  azimuths  unless  some  connection  that  ties 

L/v  and  i is  assumed, 
e 

The  displacement  seismograms  reveal  that  both  the  rise  time  and 
rupture  velocity  can  be  traded  off  to  produce  very  similar 
wave  shapes  emphasizing  that  it  would  be  difficult  to  separate  the  effects 
of  rise  time  and  rupture  velocity  using  displacement  seismograms.  Any  dif- 
ferences in  displacement  seismograms  for  various  rupture  parameters  appear 
to  be  subtle  except  for  an  infinite  rupture  velocity  where  differences  in 
rise  time  produce  noticeable  differences  in  the  pulse  shape,  as  shown  in 
Figure  8.  We  therefore  examined  the  vertical  component  of  particle  ve- 
locities and  accelerations  at  stations  B and  C to  evaluate  their  appro- 
priateness for  learning  about  details  of  the  rupture  process.  One  intuitively 
expects  that  more  i^-sight  into  the  details  of  the  faulting  are  carried  in 

the  high  frequencies  which  are  essentially  filtered  from  displacement  seis- 
mograms . 

Figvire  9 sliows  the  vertical  component  of  the  particle  velocities  gener- 
ated at  station  B for  various  combinations  of  rise  time  and  rupture  velocitv. 
Here  it  can  be  noted  that  for  a source  rise  time  of  0.5  seconds  increasing 
the  rupture  velocity  from  2 to  3 km/sec  produces  significant  changes  in  the 
amplitude  and  wave  shape  of  the  particle  velocity  pulse.  In  particular,  the 
amplitude  increases  and  tlie  wave  period  (here  taken  as  the  peak  to  peak  in- 
terval between  wave  maxima)  decreases  as  the  rupture  velocitv  approaches  that 
of  the  siiear  velocity  (3.2  km/sec)  of  tlie  medium. 
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For  a slower  rise  time  of  1.0  seconds  ttie  particle  velocity  pulses  are 


somewhat  simpler  in  character  although  the  same  gross  effect  can  he  seen. 

The  comparison  of  panicle  velocities  computed  with  and  without  the  presence 
of  a free  surface  reveal  that  even  though  the  gross  character  of  the  velocity 
pulse  appears  similar  one  can  note  differences  especially  in  the  amplitude 
of  the  first  positive  upswing  relative  to  the  second  upswing.  However,  if  one 
compares  peak  to  peak  maximum  particle  velocities  inclusion  of  a free  surfnce 
produces  an  increase  ranging  from  a factor  of  2 to  2.3. 

Figures  10  and  11  show  the  particle  accelerations  at  station  B for  three 
different  rupture  velocities  and  source  rise  times  of  0.5  and  1.0  seconds. 

The  synthetic  accelerograms  exhibit  much  character  and  individual  contributions 
to  the  acceleration  signals  such  as  the  S-wave  stopping  phase  (Savage,  IdhSi 
can  be  seen.  For  a source  rise  time  of  0.5  seconds  we  have  indicated  the 
expected  arrival  time  of  the  S-phase  and  the  stopping  phase  of  the  S-wave. 

A small  P-wave  stopping  phase  can  be  recognized  when  the  rupture  velocity  is 
3 km/sec  but  cannot  be  noted  for  slower  rupture  velocities  because  its  time 
of  arrival  occurs  after  the  onset  of  the  S-phase. 

The  accelerograms  are  much  simpler  in  character  when  the  rise  time  is 
increased  to  1.0  second  (Figure  11).  In  this  case  individual  stonping  phases 
are  not  obvious  inasmuch  as  the  rise  time  is  comparable  to  the  fault  length 
divided  by  the  rupture  velocity,  i.e.,  the  fault  is  still  radiating  energy 
from  near  the  initial  point  of  rupture  as  the  rupture  reaches  the  end  of  the 
fault.  Figure  12  shows  accelerograms  computed  at  station  C located  in  the 
backward  quadrant  to  the  direction  of  rupture.  In  these  acc-elcrograms  di  'net 

P and  S wave  stopping  phases  can  be  noted. 

The  accelerograms  presented  here  indeed  emphasize  that  the  adopted  ris# 
time  of  the  source  and  the  rupture  velocity  do  plav  an  impi'rtant  part  in  ci’n- 
trolling  the  amplitude  and  overall  shape  of  computed  near-field  acci' 1 e ri'gr  am- . 
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WhoCher  svntholic  accelerograms  are  justified  in  practice  for  comparison  with 

observed  accelerograms,  liowever,  would  seem  to  depend  on  the  azimuthal  coverage 

of  observations  availa()le,  an  assessment  of  the  importance  of  layering  and 

atteiiuation  and  whether  enough  fault  parameters  are  independent Iv  known  so  that 

an  appropriate  kinematic  model  can  be  defined. 

The  near-field  seismograms  presented  here  have  r/h  ratios  much  less  than 

5 and  there  is  no  evidence  of  a Rayleigh  wave.  Pekeris  and  Li f son  ( 1917)  have  shown 

for  a buried  point  source  that  the  Rayleigh  pulse  begins  to  emerge  at  r/h  = 

5 and  is  clearly  recognizable  at  r/h  = 10.  To  examine  the  applicability  of 

this  conclusion  for  a buried  strike-slip  source  we  computed  the  vertical  and 

radial  displacements  at  increasing  epicentral  distances  (Figures  13  and  14). 

In  these  calculations  a point  source  and  a source  rise  time  of  0.5  seconds  were 

assumed  together  with  a source  potency  U dS  of  1 x 10  unit  faults.  (One  unit 

Cl 

2 

fault  = 1000  m - km  where  dS  is  the  fault  area).  As  can  be  seen,  a Rayleigh 

* 

pulse  is  noticeable  at  r/h  = 5.  Similar  results  were  obtained  for  a dip-slip 
source  (Figure  15).  Figure  16  shows  the  effect  of  decreasing  the  depth  of  focus 
h for  various  distances  r.  Decreasing  the  depth  of  the  source  emphasizes  the 
higher  frequency  components  of  the  Rayleigh  wave  and  a large  Rayleigh  pulse  is 
prominent  on  the  seismogram  when  r/h  > 50.  The  conclusions  of  Pekeris  and 
Lifson  (1957),  and  amplified  upon  by  Mooney  (1974,  1976),  thus  also  appear  appli- 
cable to  buried  dislocation  sources. 

As  .1  lin.il  I'xample  we  examine  observed  near-field  seismograms  for  a 

= 4.t>  e.irt  h<iuaki'  recorded  at  SAGO-Kasl  and  S.Af.O-Central  as  presented  bv 
Toluison  .ind  Mct.'villy  ( 1974).  in  Figuri'  17  the  signals  shown  repri'sent  actual 
recordings  oi  ground  displacement  recorded  at  i-pici-ntral  distances  of  2.1 

* 

Results  for  a dip-slip  source  can  be  derived  using  Table  2 in  Ben-Menahcm  and 


Vered  (1973). 
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km  and  5.5  km.  We  do  not  attempt  a detailed  study  of  tlie  souree  iiaraneti-rs 
of  this  earthquake  but  rather  only  demonstrate  a fit  of  amplitmies  and  pulse 
widths  based  on  several  assumptions.  A detailed  study  of  this  earthquake 
from  a dynamic  point  of  view  is  given  by  Israel  and  Vered  (1977).  First, 
the  epicentral  location  was  shifted  by  0.75  km  in  the  directii'n  S 32°  W In 
order  to  achieve  agreement  with  the  observed  amplitudes  at  all  stations. 

An  adjustment  of  this  amount  is  within  the  uncertainty  of  the  epicenter 
location.  Secondly,  a fault  area  of  3 kra^  is  assumed  and  rupture  is  taken 
to  initiate  at  the  focus  and  propagate  to  the  southeast.  A kinematic  model 
is  taken  with  a constant  fault  dislocation  and  a constant  rupture  velocity 
of  3 km/sec.  The  particle  displacement  is  given  bv  (21)  and  the  rise  tine 
depends  on  the  distance  from  the  start  of  rupture  according  to 


where  L is  the  fault  length,  here  taken  to  be  3 km.  The  choice  of  fault 
kinematics  is  quite  arbitrary  and  the  theoretical  seismograms  shown  in 
Figure  17,  although  quite  satisfactory,  are  not  unique. 

A better  azimuthal  station  distribution  together  with  a model  which  in- 
corporates the  geological  complexity  of  the  area  could  possibly  reduce  the 
possible  solutions.  However,  details  of  the  rupture  process  are  not 
resolvable  from  examination  of  only  near  field  displacements. 

As  pointed  out  bv  McKvillv  and  .lohnson  (1974)  t hi'  comparison  o!  syn- 
thetic seismograms  with  .ictual  seismograms  provides  a method  ti'r  I'Stinatini’ 
the  source  moment  in  the  time  dom.ii  n . The  average  source  poteticv  detomiueii 
by  comparison  of  the  amplitiule  of  the  S-wave  pulse  uilh  obse r\M t 1 (>ns  at 
.SA(kl-Cen  tra  1 and  SA(k)-i;ast  is  185  cm-km'.  Convi'rting  ti'  seismic  mi'mi  nt  . 
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assuming  a rigidity  of  2.73  x 10  dynes/cm  , yields  a value  of  5 x 10 

2 

dyne-cm.  With  an  assumed  fault  area  of  3 km  this  leads  to  an  average 

23 

dislocation  of  61.7  cm.  The  computed  seismic  moment  of  5 x 10  dyne-cra 
agrees  well  with  the  value  determined  by  Johnson  and  McEvllly  (1974)  assuming 
a point  dislocation  and  the  values  obtained  by  Israel  and  Vered  (1977)  for 
a dynamical  model  of  rupture,  and  points  to  the  difficulty  in  extracting 
details  of  the  rupture  process  from  only  examining  near-field  displacement 
data  (Anderson  and  Richards,  1975).  Near  field  synthetic  accelerograms, 
however,  hold  some  promise  for  the  study  of  fault  rupture  parameters. 

An  attempt  has  also  been  made  to  examine  what  produces  the  coda  in 
the  seismic  near  field.  Theoretical  seismograms  have  been  computed  for  a 
point  source  imbedded  in  a layered  medium.  An  example  of  such  a calculation 
is  shown  in  Figure  18.  Displayed  are  the  seismograms  for  a vertical  point 
source  placed  in  a homogeneous  half-space  and  in  a layered  half-space  with 
parameters  appropriate  for  the  crustal  structure  near  the  San  Andreas  fault. 
The  source  depth  was  taken  to  be  12.5  km  and  the  epicentral  distance  5 km. 

Our  computed  seismograms  demonstrate  the  general  computational  result  that 
neither  pulse  amplitudes,  nor  pulse  durations,  are  sensitive  to  crustal 
layering  but  that  a portion  of  the  near-field  coda  can  be  attributed  to  the 
presence  of  crustal  layers  and  hence  the  receiver  structure  must  be  accu- 
rately known  before  the  source  radiation  can  be  completely  known. 
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CONCLUSIONS 


Theoretical  seismograms  have  been  presented  for  a buried,  propagating 
strike-slip  fault  to  examine  the  effect  of  various  rupture  parameters  and 
the  presence  of  the  free  surface  on  the  resulting  particle  displacement, 
velocity  and  acceleration  at  various  azimuthal  near-field  locations.  Rupture 
velocity  and  source  rise  time  usually  cannot  be  independently  assessed  on 
displacement  seismograms  inasmuch  as  the  near-field  pulse  width  is  a func- 
tion of  fault  length,  rupture  velocity,  effective  source  rise  time  and  the 
distances  from  the  beginning  and  end  of  the  rupture  plane  to  the  observing 
station.  However,  given  a judicious  (or  fortuitous)  azimuthal  choice 
of  observation  points  relative  to  a propagating  strike-slip  fault  near- 
field synthetic  accelerograms  suggest  some  promise  for  the  evaluation  of 
fault  rupture  parameters.  This  could  be  of  some  importance  in  studying 
ground  motions  in  the  near-field  of  potentially  damaging  earthquakes. 

For  the  strike-slip  model  discussed  here,  allowance  for  the  presence 
of  a free  surface  by  doubling  of  the  displacements,  velocities  and  acceler- 
ations computed  in  an  infinite  space  appears  to  be  valid  to  within  25%  or 
so.  This  conclusion  is  not  generally  valid  in  the  case  of  the  residual 
static  displacements.  A portion  of  the  near-field  coda  can  be  directly 
attributed  to  the  presence  of  crustal  layers  which  must  be  accurately  known 
before  precise  estimates  of  source  parameters  can  be  made. 
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figure  16 


FIGURE  17 


APPENDIX  A 


SAMPLE  COMPUTER  LISTING 

There  are  two  programs  to  compute  the  vertical  component 
of  motion  for  a propagating  strike-slip  fault.  Program  one 
should  be  submitted  first,  the  numerical  results  of  this  pro- 
gram being  stored  on  disc.  This  information  is  used  when 
program  two  is  submitted.  The  programs  are  used  to  compute 
ground  displacements,  velocities  or  accelerations.  The  source 
may  bo  either  a line  source  or  a rectangular  source.  Compu- 
tations may  be  done  for  an  infinite  or  half  space  model. 

The  most  time  consuming  computations  are  the  functions 
u , u , (see  eqns.  (15),  (16).  Program  one  computes 

Z,p^  Z,Sj^ 

the  numerical  values  of  these  functions  once  and  for  all.  The 

values  of  u + u are  stored  on  disc  for  all  subrec- 

z,p^  z,s^ 

tangular  sources  required.  Program  two  uses  these  values  sever- 
al times  to  compute  the  acceleration,  velocity,  or  displacement, 
for  any  required  rise  time  or  rupture  velocity. 
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PROGRAM  ONE 


INPUT  (see  figure  1) 

ALFA:  P velocity  (km/sec) 

BETA:  S velocity  (km/sec) 

H:  Source  depth  in  km  (depth  to  top  of  fault  plane) 

D:  Horizontal  range  in  km 

DT:  time  increment  in  u + u computations. 

z,p^  z,s^ 

Time:  the  required  computational  time  of  ground  motion  (sec) 

E^^:  the  absolute  permitted  error  in  u , u 

z,p  ' z,s- 

+ + 

EE  : the  relative  permitted  error  in  (same  as  above) 

DELH:  fault  width  (km) 

DELX:  fault  length  (km) 

PHI:  azimuth  to  station  in  degrees  (see  figure  I) 

NH:  number  of  subrectangles  required 

INDXl:  1 denotes  a half  space  model;  2 an  infinite  space  model 
INDX2:  1 denotes  a rectangular  source;  2 a line  source  at  depth 
H,  length  DELX. 

OUTPUT 

1/  Arrivals  time  of  P,  SP,  S from  depth  H and  H + DELH  for  all 
subrectangles . 

2/  The  numerical  values  of  u + u for  all  rectangles. 

z,p^  z,s^ 

These  values  are  printed  and  stored  on  the  disc  in  unit  8. 


++ 


The  error  less  than  the  absolute  error  or  the  relative  error. 


i 
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PROGRAM  TWO 


INPUT 

NV : The  number  of  different  rupture  velocity  cases  to  be  com- 

puted . 

INDEX:  1 for  displacement  ground  motion,  2 for  velocity  ground 

motion,  and  3 for  acceleration. 

TZ:  rise  time  in  seconds  (2  in  equation  (21)) 

VZ : the  first  value  in  rupture  velocity  (km/sec) 

DV : increment  in  rupture  velocity 


OUTPUT 

A/  The  values  of  rupture  velocity  and  rise  time  are  printed  as 

well  as  the  values  of  ground  motion.  Displacement  amplitude 

is  given  in  terms  of  U^,  velocity  in  terms  of  U^/(sec),  and 

2 

acceleration  in  terms  of  UQ/(sec  ). 

B/  A plot  for  each  required  rupture  velocity.  The  values  of 

displacements,  velocities  and  accelerations  were  multiplied 
by  1000  before  plotting. 
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PROr.RAM  ONE 


C PROGRA^<  TC  CUMPUTE  SURFACE  MOT  UN  OF  mAlFSPACE  t DuE  TO  dORlEO 

C SOURCE  VERTICAL  MOTION 

implicit  RtAL*a  (A-H,Q-£I 

DIMENSION  U(  1008  ) .US(  1 008) .UU(  1008)  .UK(  1008  ) 

COMMON/PLO/Q. US .DT.TZ 

CCMMCN/PPl /ALFA. BETA.  T IME .E.EE. I PLOT 

COMMON/AA A/INDX 1 , INDX2 

READ  70 .ALFA#  BET A . H, D . OT . T I ME. E . EE 

read  70  . OELH.OELX.PHI 

READ  73, NH 

READ  73 . 1 AOXl  , INDX2 

PRINT  104, ALFA. BETA, H.O.OT. TIME. E.EE 
PRINT  74, CELM.DELX.PMI 
PHINT75.NH 
PRINT  75  , INOXl , INOX2 
70  FORMATI8F10.5) 

73  F0RMAT(3I£) 

DO  207  I = l ,1LC0 
207  Q( I )=0. 

H = H«-0ELH/2  .DO 

WRI TE(8)  H. D. DT. TIME. DELH, del X. PHI 
WRITEta)  NF, ALFA, BETA 

74  FORMAT! IX . 5El 5.7 ) 

PI=DARCOS  <-l  .DO ) 

PHI=PHI*P I/l .8D2 
D 1 =D 

RNS=D14Q1 

RS=RNSF(H-DElH/2  ,00)*(  H-DELH/2  .0  0 ) 

R=DSURT(RS  ) 

T1=R/ALFA 
IP  = Tl/OT-t-  1 
NI=T1ME/DT  VIP 

75  FORMAT! IX .513 ) 

H1=M 

ANM=NH 

OXT=l,O0/ANH 
ALC=0.O0 
TS=0  . 

DO  27  11  = 1 .NH 
AI I=II 
TS=TSfOXT 
AL=DELX*! TS) 

ST=AL-ALC 

SP  = ! ALVALC  )/2 .DO 

ALC=AL 

RNS=01*Dl  «-SP*SP-2.DO*Ol4SP*DCOS!  PHI  ) 

DIM=OSORT!RNS) 

H=H1  ♦!  .5D0*Dc.LH 

IF! I N0X2.EQ.2 )H =Hl F 0 . 5 O0*DEL H 
D = OIM 
MNH  = 2 

IF ! I N0X2. e U.2 )MNH  = I 
DO  2 1=1. 1C08 

2 QQ! 1 ) = 0.00 
DO  I J=1.MNH 
A J=-  I . 

1F!INDX2.EC.2)AJ=1. 

IF! U.EU.2  ) AU=I. 

H=H-0ELH 

301  call  VRM!H.0.M) 

DO  I i=t.icoa 
QQ!  I ) = QQ1  I )FU!  I M>AJ 
1 CONTINUE 

DO  3 1=1.1 008 

3 O! I ) =UU!  1 ) 

WRITE (8  ) ! Q(  1 ). I =1 .NI  ) 

print  104  . (Ui  I ) . I =1  .NI  ) 

27  CONTINUE 
104  format !2X . lOE 12 .4 ) 

STOP 

END 
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subroutine  VR*i(  HtOlH.NI  ) 
implicit  RtALOa  (A-H.O-Z) 

EXTERNAL  F UNP  .FUNS  .FUNSP  .fJNSS 

DIMENSION  UUOOat  . U S ( 10  0 8 > . USK  1 00  8 ) . 0 S2(  1 0 OB  I 

COMMON/PLO  / U.QS.OT.TZ 

COMMON  HR2 .0R2. SA* SB, SC, RA .RB, TBP . TB S.U  ,T  . TbSR 

COMMQN/PPl/ALFA.OETA.  TIME, £. EE, IPLQT 

COMMON/ A A A/INOX  1 . IN0X2 


,US3( 1008) 


0=01  M 

00  999  1 = 1 .1008 
0(11=0.00 
QSl ( I )=0.00 
QS2 ( I )=0 .00 
0S3< I ) «0.C0 
999  0S< I ) = 0.00 

PI K=O.SOO*CARCOS(-1 .00 ) 

PRINT  73,H,0 

73  FOHMATOOX. ‘VERTICAL  DISPLACEMENT*./,  18X,  ‘SOURCE  DEPTH  ‘,F8.3, 
1 ‘KM  . horizontal  range  ‘,F8.3,‘  KM,  ‘ ) 

RS=  HAH-*-D*0 
0R2=0/RS 
HR2=H/RS 
R =DSORT(rtS) 

FA  =1. 00/ALFA 
EB  =1. 00/BETA 
SA  = FA*FA 
SB  = EB«£B 

HCR=H/DSORT(SB/SA-l .00 ) 

SC  = O.SASe 
RA  = RS«SA 
R8  = Rs*se 
TBP  =R*FA 
PRINT  74, TBP 

74  FORMAT! 18X  ,‘P  ARRIVAL  TIME  ‘.F10.5  ,‘  aEC,‘) 

IP  =TBP/DT  ♦! 

TB  =IP4DT 
T=  TB-OT 
N=NI 

DO  1 I=IP.N 
T=T»DT 

1 Q( 1 ) =QUAC( O.OO.P IK.FUNP.E.EE .4 ) 

TBSR  =R*EB 

TBSH  =0*FA  ♦H40SQRT(SB-SA) 

IF(lNOXi.£0.2>  GO  TO  3 
IF(O-HCR)  3.3,4 
3 TBS  =T8S« 

75  FORMAT!  18X  ,‘ NO  SP  HEAD  WAVE.  S ARRIVAL  TIME  ‘,F10.S,‘  S£C.‘) 

PRINT  75, TbSR 

IS  = TBS/DT 

J = iS-IP 


TS= IS*OT 
T = TS-DT 

00  2 I=1S,N 
T=T*  DT 

QS!I)=  QUAD! 0.00, PIK.  FuNS,E.EE,4) 

2 0! 1 I =0! I ) «QS! I > 

GO  TO  10 
4 TBS  =TBSH 

PRINT  76. TBSH, TBSH 

76  FORMAT! 18X.‘  SP  HEAD  WAVE  ARRIVAL  TIME  ‘,F1C.5,‘  SEC.  3 ARRIVAL  TI 
IME  ‘.F1C.£.‘  SEC‘) 

IS  = TBS/DT  *^1 

TS=I S*0T 
J = IS-IP  *l 

T = TS-DT 

PORH=OARSIN  (OSORT! T8SH/TBSR>  ) 

DO  18  I = IS,N 
T=T»0T 

IF!T .GT.TBSR)  GO  TO  9 
PCRH=0ARS1N  ! OSQRT! T8SH/T » ) 

QSl!  I ) = QUAO!PLRH  ,PI K,F  UNSP ,E ,EE, 4 ) 

Q!  i > =Q! 1 ) ♦USl  ! I > 

GO  TO  18 

9 QS2( I )=QUA0!0 .00 ,PI K . F UNS . E , EE .4 ) 
QS3!I)=auAD(P0RH.PIK,FUNSS.E.E£,4) 

Q!  I ) =Q! I ) ♦053! 1 ) ♦QSZ!  I ) 

18  CONTINUE 
10  FACT0R=1. 

RETURN 

END 
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FUNCTION  FUNPIX) 
implicit  AcM-*a  lA-H.O-Z) 

COMPLEX*! e C . V . VR. VSQ.P.TVP .0 .VR2. AV2. bV2 .CDSQWT 
COMMON  liR2.0R2.SA,  SB. SC. RA  .RB  . TBP  .TBS  , RANGE  . TT 

COMMON/ AAA/INOXl . IN0X2 
C=:(0.00.1  .00) 

T=TT-(TT-TaP  ) *DSIN(  X ) * *2 
C IN  THIS  SUB  T MEANS  TAU 

A=  T*T-RA 

IF(A.LT.O.OO)  A=0.00 
A=  OSQRTIA) 

V=DR2*T*HH2*A*C 
VSQ=  V *V 
P =SC  - VSU 
AV2aSA-VSQ 
BV2-SB-VSC 
VR=V*RANGE 
DT=TT-T 
TVR  = DT*WR*  VH 
DT=CT*OT 
VR2=VR* VR 
0=CDSQRT( AV2) 

C=0 

I F( I NOXl .E0.2  ) GO  TO  1 

0=V*P*0/( ( P*P*VSQ*D*C0SURT<aV2 ) ) *COSQRT ( TVR *( T*TBP) I ) ♦( DT*T VR* VR2 ) 
GO  TO  2 

1 0=V*0  /(SB  *CDSQRT(TVR*(T*T8P)I )*(0T*TVH*VR2 ) 

2 1F( lNOX2.EQ.! ) 0=D/C 
FUN*0 

FUNP  =FUN 

RETURN 

END 


FUNCTION  FUNS(X) 

IMPLICIT  REAL*8  (A-H,Q-2) 

COMPLEX*! 6 C .V  . VR.VSO.P.TVR .D .VR2. AV2. BV2.CDSQRT 
COMMON  HR2 ,DR2.SA.Sa.SC*RA.RB.TBP .T BS • RANGE. TT . T BSR 

COMMON/ AAA/INOX ! . INOX2 
C=( 0.00 . ! .00) 

T=TT-(TT-TBSR )*OSIN<  X) **2 
A=  T*T-RB 
IF (A.LT.O.  ) A=0.00 
A=  OSORTIA) 

C IN  THIS  SUB  T MEANS  TAU 

V=DR2*T*HR2*A *C 
VSO=  V *v 
P a sc  - VSU 
AV2=SA-WS0 
BV2=Sa-VSQ 
VRaV*RANGE 
OTaTT-T 
TVR=OT*VR* VR 
DT=OT*OT 
VR2= V«*WR 
0=COSQRT( AV2) 

CaCOSORTI eV2) 

IF(INOXI.E0.2)  GO  TO  ! 

Oa-V*8V2*(OT*TVR*VR2 ) *0/( (P*P*  VSQ*0*C0SURT ( BV2) > *COSORT ( T V R* ( T * 
*TBSR  ) ) ) 

GO  TO  2 

! 0=-V*C  * ( 0T*TVR*VR2 ) /<SB  *COSOR T( T VR* ( Tt 

*TBSR)  ) ) 

2 IF( IN0X2.E0.!  ) QaO/C 
FUNaD 
FUNS  =FUN 
RETURN 
END 
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FUNCTION  FUNSP(X) 

IMPUiCIT  HtAC*8  (A-H,Q-^^ 

COMPLEX«ie  AV2.0.C0SQRT 

COMMON  HR2.0R2*SA.Sd*SCt  RA .R3 .TOP  «TBS .RANGE#  TT 

COMMUN/AAA/lNCXl . 1NOX2 
T=TT  *OSIN(X)**2 

A=  T»T-Rd 
A=OSaRT(-A) 

C IN  THIS  SOU  T MEANS  TAU 

V=DR2*T-HR2*A 
WSU»  V *V 
P =SC  - VSO 
AV2=SA- VSC 
BV2=Sd-WSQ 
VR=V«RANG£ 

OT=TT-T 
TVRsOT^WR  + VR 
OT=CT*OT 
VR2=VR*VR 
U=COSURT( AV2) 

D=V*aV2*(OT*TWR«-VR2)  *D»  OSQRTI  T/  ( T VR*  < T «-T  BSR  ) *(  TBSR- 
«0«COSaRT(BV2«(l .00*0.001  1 I 
IF ( 1N0X2.EQ.1  1 D=0/C0SQRT(BV2«  < 1 .00 .0.0  0)  ) 

FUN^OVI 0.00.- 1.00) 

FUNSP=  FUN 
RETURN 
END 


FUNCTION  FUNSSIX) 
implicit  real *8  (A-H.O-Z) 

C0MPLEX«16  AV2.0.COSQRT 

COMMON  HR2.0R2.SA.SB.SC.RA.R6.TBP .TBS .RANGE. TT.TBSR 

COMMON/AA A/INOX  1 .IN0X2 
T=TBSR  *0SIN(X)**2 

A=  T»T-RB 
A=0S0RT(-A) 

C IN  THIS  SUB  T MEANS  TAU 

V=DR2*T-HR2*A 
VSO=  V *w 
P =SC  - VSQ 
AV2-SA-VSQ 
8V2=SB-VSQ 
VR=V«RANGE 
OT=TT-T 
TVR=OT>VR  ♦ VR 
DT=OT*-OT 
VR2= VR*VR 
D = CDSQRTt  AV2) 

0*V*aV2*IDT»TVR*VR2)*0*  OSQRTI T/ < T VR* ( T ♦TBSR 1 • ( T T -T) ) )/( P*P*VSQ* 
*0*C0SORT( e V2*( 1 .00.0.001)1 
IF( IN0X2.EU. 1 1 U-O/COSQRTI BV 2* ( 1 . DO  * 0 .0 0 1 ) 

FUN=0*(0.00.- 1.00) 

FUNSS=  FUN 
RETURN 
ENO 


.TBSR 


T) ) )/IP*P*VSU* 
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FUNCTION  aUAU(A.B.<^UN.£PS.ETA.MI  N) 
implicit  RtAL«e  IA-H,0-Z) 

OIMENSICN  Q(16) 

FC^^=FUNI  A) 

FCNd  =FUN<8» 

TABS  =OA6S<H  ) *(OABS(FCNA)«-OAOS  (FCN9)  ) *0.5 

T = H •(  FCNA  ♦ FCNB  ) *0.5 

NX  =I 

00  12  N=1 • IS 

H=H*0.5 

SUMsO.OO 

SCOAA^O.OO 

SUMA8S=0.00 

00  2 I^l.NX 

XI-2.*  OFLOArtll-l. 

FCNXI  = FLN(A*XI*H) 

SUMABS>SUMABS*  OABSI FCNXI) 

FCNXI  = FCNXI  ♦ SCORR 
SS  = SUM  *FCNXI 
SCORR-=  SUF-SS  «■  FCNXI 

2 SUM=SS 

T=  T*0.5*H*SUM 

TA8S=  TABS«O.S  * CABS(  H ) «SUMABS 
Q(N)=  2.* ( T*H*SUM)*0. 33333333 
IF(N-2)  10.3.3 

3 F =4 . 

00  4 J-2.N 
I =N* I- J 
FsF*4. 

4 Q(1)=Q(1«-1)*(0(I*1>-0(I))/(F-1.) 

IF(N-3»  9.5,5 

5 IF(N-MIN)  9,6.6 

6 K-  DABS(U(  l)-QX2)4-OA8S(UX2-QXl  ) 

IF  (TABS)  7,8.7 

7 IF(X/TAaS-3.*(OA8S<eTA)*0. 149011610-7))  11.11.8 

8 IF(  X-3.«0AbS(EPS)  ) 11.11,9 

9 QX1=QX2 

10  QX230( 1 ) 

12  NX=NX*NX 

1 1 QUAD  - Q( 1 ) 

RETURN 

END 
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fOO 

200 


Sl.uKL  Jl  INl  oKA^Ji  (N.hAC  , am  I , YMI  ,aMA,  YMA) 

Cl  'li-lL.i'i/ bi"  / lOOojfAPt  IJJJ)  f uJI-l-  1200) 
CU(-il'1Ci</Pl'  / Awt  iJOd) 

UT  =r AC 

T=U. 

Di)  ^00 
XP i 1 ) =T 
T=  T4-UT 

UU  2o0  J=i,N 

YH  i J )=A^l J )*1CU0. 

PC^^^d'♦)  = l • 
r)Ut"("  I i 02  ) = 1 . 
lUjt-rll/7)=i. 

CA»-L  CoJcCcluoFFtA'II  fY*^!  fXMAiYi'IA  i 
call  ^uujCijIOoFF  jO^t^OO^fd#!  ’j0») 

CAL*.  vjKAPFLldLirF»NiXFfYr’»i»'A'  flt'2'tOf  *X 
PE  T OK iM 

E N u 


V S2'  ) 


SCuKLufii^c  bl  iT2»ur»  *) 

RfAL*d  L bf  jTiJb  , TAOd 
CUI  Mln  /KlL/  l»  CSf  jTJlj  » E^-ub 
C -Km'IOi'^/  i N / 1 H Ct  A 
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